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Abstract. The use of the QR factorization of the Hankel matrix in solving the partial 
realization problem is analyzed in this paper. 

Straightforward use of the QR factorization results in a new realization scheme that 
possess all of the computational advantages of Rissanen’s realization scheme. These latter 
properties are computational efficiency, recursiveness, use of limited computer memory, and 
the realization of a system triplet having a condensed structure. Moreover, this new scheme 
is robust when the order of the system corresponds to the rank of the Hankel matrix. 

When this latter condition is violated, an “approximate” realization could be de- 
termined via the QR factorization. In this second scheme, the given Hankel matrix is 
approximated by a low-rank non-Hankel matrix. Furthermore, it is demonstrated that 
column pivoting might be incorporated in this second scheme. 

The results presented in this paper are derived for a single-input/single-output sys- 
tem, but this does seem not to be a restriction. 

Keywords: Partial realization, QR factorization, Hankel matrix, Markov parame- 
ters 
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1. Introduction. The so-called Hankel approach [l] to solve the partial 
realization problem [2] is studied in this paper. The cornerstone of this approach consists of 
the following two properties of the Hankel matrix, constructed from the Markov parameters: 
(1) the so-called shift invariance property and (2) the correspondence of the ‘'numerical” 
rank of the Hankel matrix with the order of the system. Basically from these two properties 
“any” factorization of the Hankel matrix allows to solve the partial realization problem. 

This observation was first made by J. Rissanen [3] in 1971. The algorithm described 
in [3] possesses a number of appealing properties: it is recursive and, computationally as 
well as storage related, very efficient. However, the original presented scheme was not 
numerically stable. This was observed by L.S. De Jong [4], who in the same paper proposed 
algorithmic modifications to make the scheme numerically stable. 

The remaining drawback of this computational scheme was that it only used a very 
small amount of the available Markov parameters. Hence, making it only applicable when 
the Markov parameters are infinite-accurate. For the latter circumstances, which of course 
are very unrealistic, Rissanen’s realization scheme will be indicated in this paper as an 
elegant solution, summarizing the above mentioned algorithmic properties. Additionally, 
this realization scheme directly resulted in a state space description in a canonical form [3], 
[4], what will be indicated as an attractive property of the realization scheme. 

To make the Hankel approach applicable for practical realization problems, S.Y. 
Kung proposed to use the singular value decomposition (SVD) to factorize the Hankel 
matrix [l]. This allowed first to reliably determine the rank of the Hankel matrix and 
second to consider all available Markov parameters. Although the SVD might lead to a 
numerical robust solution of the partial realization problem, it is still, for large systems 
and/or marginally stable systems, very storage and computational inefficient. Also, it does 
not result in a special canonical form and therefore certainly does not possess the elegantness 
and attractiveness of Rissanen’s numerical stable scheme. 

An inherent drawback in the Hankel approach is that even (relatively) small errors 
on the Markov parameters cause the second basic property to be violated. This is due to the 
repetition of the individual Markov parameters and hence, so do their corresponding errors. 
An illustration of this phenomena will be given in this paper. Therefore, only factorizing 
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the Hankel matrix to solve the partial realization problem is not sufficient, no matter which 
factorization technique is used. This drawback has stimulated additional research which 
has been indicated as reduced-order modelling. Here one focuses on “approximating” the 
given Hankel matrix by a low-rank Hankel matrix [5], [6], and [7]. For certain practical 
applications [8] the latter constraint might be relaxed to approximating the given Hankel 
matrix by a low-rank non-Hankel matrix. Also here the key tool used in solving this type 
of realization problems is the SVD. 

The incentive of this paper is to investigate the usefulness of the QR factorization of 
the Hankel matrix in solving partial realization problems. Two types of realization problems 
are investigated: (1) when the two fundamental properties of the Hankel matrix, mentioned 
above, hold, and (2) when the given Hankel matrix is approximated by a low-rank non- 
Hankel matrix. 

The organization of this paper is as follows. After (re)stating the partial realization 
problem in section 2, the QR factorization is first analyzed for the type of partial realization 
problem, where the two basic properties are not violated. In this case, the Markov param- 
eters are indicated to be “accurately” given. Section 4 then discusses the use of the QR 
factorization in solving the realization problem from a low-rank non-Hankel approximation 
of the given Hankel matrix. This is indicated as “approximate” reduced-order modelling. 
The presented two realization schemes in sections 3 and 4 are then evaluated experimentally 
via a number of numerical examples in section 5. Finally, section 6 present the concluding 
remarks. 

The discussion throughout this paper is restricted to the single-input/single-output 
(SISO) case. But the results rather straightforwardly extend to the multi-input /multi- 
output (MIMO). 

2. The partial realization problem, in this section we outline 
the partial realization problem and indicate the two fundamental properties of the Hankel 
matrix, which form the cornerstone in solving this problem. This was also indicated in [l] . 

Let be the impulse response of a linear time-invariant system: 

**+i = Ax k + bu k (1) 
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Vk = c T x k (2) 

where x * 6 i? n . Generally, the terms /i t are called the Markov parameters . Given the system 
description (1-2), they could be stated as follows: 

hi = c T A i ~ 1 b (3) 

The so-called partial realization problem is one of recovering the system triplet {A,b,c T } 
from the finite sequence {/i,} for t = 1,2,* •• M + JV - 1. Here the integers M, N have to 
be larger than n, however since n is not known a priori they are generally (assumed to be) 
taken much larger than n. 


The Hankel matrix H q (k,i) that is used to solve this problem is constructed from 
from the Markov parameters as follows: 


^ hq hq + 1 ^q+t — 1 


H q (k,e) = 


h 


?+l 


^ hq+k - 1 


hq+k+e - 2 J 


( 4 ) 


The two fundamental properties of this Hankel matrix are: 


Property 2.1: The rank of the Hankel matrix Hi(M,N), taken M and N “sufficiently 
large” as indicated above, is the order of the state space system given in (1-2). 


This property can easily be understood by using the Cayley-Hamilton theorem [2] and 
factorizing H\(M,N) as: 


/ 





\ c t A m ~ 1 



( 5 ) 


where 0 6 R Mxn and C € R nxN . 


Property 2.2: Using a similar factorization for the Hankel matrix the shift- 

invariance property of this Hankel matrix is: 

H 2 {M,N) = OAC (6) 
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Remark 2.1: From (5-6) we directly observe that any factorization of the Hankel matrix 
as given in (5), allows to determine a state space realization. The advantage of using the 
SVD is twofold: first, the calculation of the system state matrix A , via (6), can reliably and 
easily be done once the SVD of the Hankel matrix is available and secondly, the order of 
the system n is also reliably revealed by a SVD. 

However even for (relatively) small errors on the Markov parameters, the information 
about the order of the system, which is crucial in the solution of the partial realization 
problem, cannot be retrieved from the rank of the Hankel matrix. For these circumstances, 
which will be defined precisely in the sequel, properties (5-6) are no longer sufficient to solve 
the partial realization problem. Therefore, additional conditions are required [4], [5]. 

In the next section we start treating the realization problem where both fundamental 
properties hold. 


3. Accurate Markov parameters. 


3.1. The order of the system given. Let us first assume in this section that 
the Markov parameters are exactly known and that the order of the system (1-2) n is also 
given. Then, consider the QR factorization of the Hankel matrix Hi(M,n + 1): 


H 1 (M,n + l) = (Q 1 |Qs) 




V u J 


QiR 


(7) 


Here M > n and Q\ € R Mxn with Q\Q\ = In- Furthermore let us introduce the following 
notation: 


Q 1 € R( M ~ 1 ) xn — denotes Q\ without its last row 
<5i — denotes Q i without its first row 

/? g i? nx ( n+1 ) — is an upper trapezoidal matrix given as [rj • ■ ■ r n+ i] 

The fact that R is upper trapezoidal and of rank n will be shown later on in corollary 3.3. 


Using the above notation, the solution to the partial realization problem becomes: 
cj = (g n </i 2 ■ ■ * 9in) denoting the first row of Qi (8) 

b c = ri (9) 
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and A c can be computed from either: 


Q 1 A C = Q i (10) 

or, 

Ae(ri--T B ] = [r 8 --T n+1 ] (11) 

The solution, given by equations (8, 9 and 10 or 11) is completely similar with the one base 
on the SVD [1]. However, now the algorithmic complexity is less. For example, when using 
Eq. (10) to compute A c , Q i only has to be accumulated and this does not require iteration 
as in the SVD to compute the left singular vectors. 

In the decomposition (7) [rj • • ■ r n ] can be considered as the controllability matrix of 
the realizing pair { A e , b c }. This then should imply that this pair is in controller-Hessenberg 
form. This assertion is proved in the following theorem. Let us first define the controller- 
Hessenberg form [9]. 

Definition 3.1: When the compound matrix ( B | A) of the system triplet {A, B,C} is 
upper trapezoidal, the state space description {A, B,C} is in controller-Hessenberg form. 

Theorem 3.1: The system triplet given by Eqs. (8, 9 and 10 or 11), obtained via a QR 
factorization of the Hankel matrix given in (7) is in controller-Hessenberg form. 

Proof: From definition 3.1 for SISO, the controller-Hessenberg form, requires A e to be in 
Hessenberg form. This can directly be concluded form (11) since [r^ • • • r n +i] is Hessenberg 
and [fj • • • r n ] is upper triangular. 

This in combination with the structure of the input-vector b c given in (9) defines the 
controller-Hessenberg form. 

From this theorem, we can now derive a constructive algorithm to solve the partial 
realization problem with the order of the system given without explicitly accumulating the 
used orthogonal transformations in the QR factorization of the Hankel matrix. 

Algorithm 3.1: 

1. Construction of A c , given as [aj •••<*„] from (11): 

Ol = r 2 r U 
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for i == 2:n 
end 

2. Construction of cj from: 

cH r l • • * r n] = \ h l • * * K\ (13) 

3, Construction of b c given by (9). 

Corollary 3.1: Using the RQ factorization of the Hankel matrix 1 , M) for M > n, 

results in a similar way into the observer-Hessenberg form, defined in [9]. 

Remark 3.1: Since Eq. (13) can straightforwardly be executed recursively , algorithm 3.1 
is a recursive scheme to realize a system triplet {A c , B c ,Cj} once the Hankel matrix has 
been triangularized. And the latter operation can also be done recursively. 

Remark 3.2: In [3], the observer-Hessenberg form was also derived without explicit notice. 
However, the developed implementation was numerically unstable [4] and furthermore the 
realization was obtained only based on a very restricted number of Markov parameters. 
This latter drawback may even for “very small” errors on the Markov parameters, result in 
“large errors” on the realization. This is indicated explicitly in the experimental evaluation 
study in section 5.1. These two restrictions do not apply to algorithm 3.1. 

Remark 3.3: Algorithm 3.1 also applies to the MIMO case. For this case, we recognize two 
situations. First, the dimension of the input is a multiplicity of the state dimensions. Then, 
algorithm 1 straightforwardly holds, except now that the scalars and vectors a t in (12) 
are respectively square and rectangular matrices. Secondly, when the above multiplicity 
does not hold, the algorithm remains similar, except for the computation of the last row of 
A c . This requires the use of pseudo-inverses [10]. 

3.2. Detecting the order of the system. Generally, the order n of the state 
space model is not a priori known. Based on property 2.1 it can however be retrieved from 


a, = ( **1+1 - \ai • - 


( r, ' 


V r (--i)« J 




-1 

a 


( 12 ) 
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the Hankel matrix. Particularly, the QR decomposition of this matrix directly displays the 
order of the system. In order to show this, let us first state the following theorem. 


Theorem 3.2: When the state space system (1-2) is of order n, then the n+ l,n + 2,- • • , N 
columns of with N > n are linearly dependent on the first n-columns. 


Proof: Let the characteristic polynomial of the considered n th order system (1-2) be given 
as: 

a(z) = z n + d\z n * + • • • + d n -\z 4- d n (14) 

Then according to the Cayley-Hamilton theorem [2], we can write the (n -I- l) th column of 

as: 

/ h n+1 

hn+2 


= Hi(M, n)d 


(15) 


V hn+M J 

where d = ( d n d n -i • • • di) T ± 0. This completes the proof for the (n -I- l) th column vector 
of Hi(M, N). The ( n + 2) th column of Hi(M,N) can then be written as: 


+ d\Hi(M,n)d— H\( M, N)d x (16) 


/ h \ 


0 ) 

h n +Z 


dn 


= H 1 (M,n ) 

dn— 1 

^ ^n+M+1 j 


d 2 j 


with d 1 ^ 0. 


In a similar way the n + 3,- • • ,N columns can be expressed as a linear combination 
of the first n columns of H\(M, N), what completes the proof. 

Corollary 3.2: Based on property 2.1 and theorem 3.2, the first n columns of 
are linearly independent. 


Corollary 3.3: When the Markov parameters are exactly given, the QR factorization of 
becomes: 


H 1 (M,N) = Q 


Rn | Rn 


n 

N - n 


(17) 
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with Q e r M xN . 

Hence, a QR factorization without column pivoting of a Hankel matrix of sufficiently large 
dimensions allows us to directly read-off the order of the system. 

Generally, the conditions on the Markov parameters as imposed in corollary 3.3 are 
not met in practice. Therefore, we are interested for the size of the errors on the Markov 
parameters for which Corollary 3.3 holds. This is investigated in more detail in the next 
section. 


3.3. Perturbation analysis. Let us assume that the Markov parameters /i t , 
defined in (3), are perturbed by a data acquisition error e t , which satisfies: 

hi < o (18) 

Hence, the Hankel matrix is now given as H\(M, N) + E, where the perturbation matrix E 
is a Hankel matrix itself. In this section we additionally assume that these perturbations 
are acute. This latter notion is taken from [11] and is defined as follows. 


Definition 3.2: The matrix (A + E) is an acute perturbation of A if and only if: 


rank(A) = rank(A + E) 


(19) 


This condition implies that || ([ 2 1| £7||2 < 1 [10], where A^ denotes the pseudo-inverse of A 
[11] and ||(.)||s = SU P|| X || 2 = lllOlh- 

Under the presence of the perturbations defined above, the QR factorization given 
in (17) becomes: 


( Hi(M, n) + £1 | H n+ \(M, N - n) + E 2 ) ={Q + W) 


Ru + ^1 I #12 + F2 


0 


R 22 


( 20 ) 


This decomposition clarifies that theorem 3.2 no longer holds “exactly”. For this case it 
is however more natural to refer to the linear dependency in least squares sense. Let us 
therefore define the following least squares problem associated to the decomposition (20). 


+ Ei)x - (H n+1 (M,l) + E 2 )\\ 2 (21) 

Here we focus on the way the (n + l) th column of N) + E) is lying in the space 

spanned by its first n columns. This information is incorporated in the residual of (21). 
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Therefore, the task now is to examine the influence of the perturbations {E\, E 2 ) on this 
residual. The actual influence is of course very difficult to calculate. However based on [11], 
we can formulate an upperbound for the perturbation on this residual due to the errors 
{Ei, Ei). This can be retrieved from the following lemma, taken from [11], 

Lemma 3.1: Let the residual of the following least squares problem: 

min*||(A + E\)x - (b+ 22 2 )||2 


be denoted by r and the actual residual of the unperturbed least squares problem min z ||Ai- 
6) ]|2 be denoted by r, then an upperbound for the error on r is: 

||(A+Bi)t 11,11^11, 


||r-r||, < \\Ehh + - . 

[1 + (||(A + -Ei^lhll-Eiih) 2 ] 5 

The upperbound in (22) can for the errors under consideration be taken as: 

II? - r ||2 < II-E2II2 + IIOA + £i)t|| 2 ||£i||2 


( 22 ) 


(23) 


Assuming now that the original system was of order n, the application of lemma 3.1 yields 
the following upperbound for the residual r 22 of (21): 

|r«| < (1 + V^\\(Hi(M,n) + E^\\ 2 )ay/M (24) 


Therefore, by calculating this upperbound for different values of n, we just have to inspect 
the magnitude of the diagonal elements of the R-factor in (20), relative to these bounds, in 
order to detect the order of the system. This can be incorporated during the calculation of 
the QR factorization of the Hankel matrix, a s is demonstrated by algorithm 3.2. 

Remark 3.4: The diagonal elements ra of the R-factor in (20) can be expressed into the 
following product: 

Ml) ■ a c(2, 1) • • •a c (i,t - 1) 

what exactly corresponds to the controllability measure, defined for SISO in [12]. In this 
perspective, a small r,, indicates the inclusion of a “weakly” controllable mode in the real- 
ization. The performed error analysis allows to precisely define the notion “weakly” here. 
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Algorithm 3.2: 


tf 2 ° = H\(M, N) 
for i = 1 : N 


Ii- 


Q, j 




0 I HI 


t-1 


' gj i « j 
o i «i ) 


( 25 ) 


(where /, denotes the identity matrix of order i (with /o = { }) 

H\ e R ixi 

H' 2 e /j(M-i)x(W-i)) 


if ||ffj(:,l)|| 2 <(l + (26) 

(where H 2 (:,l) denotes the first column of H 2 ) 
then go to OUT: 

end 

OUT: n = « 

Define the first (t + 1) columns of [ H\ \ ★] in (25) as [rj • • • |r* +1 ] which is then 

used as an input for algorithm 3.1. to compute the triplet {A c ,b c ,cJ} 

Remark 3.4: The quantity || (#0^112 in (26) can be computed by the inverse iteration 
method [13]. The iterations involved here are very efficient since H\ is upper triangular. 

Remark 3.5: The calculation of ||Rj(:,l )||2 in (26), each step of the do-loop in algorithm 
3.2, can be done very efficiently and recursively. This is because each two adjacent columns 
of the Hankel matrix only differ in their first and last element. 

4. Approximate reduced-order modelling. Because of the rep- 
etition of the Markov parameters and hence also their errors, the practical applicability in 
applying algorithm 3.2 straightforwardly might be limited. 
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The limitations are that the actual size of these errors on the Markov parameters have 
to be “relatively” very small in order to satisfy condition (19). Of course depending on the 
size of the dimensions of the taken Hankel matrix, the assumed upperbound a in (18) might 
still vary. However, here we are facing another intrinsic limitation of the Hankel approach. 
Namely, that when the accuracy of the Markov parameters is “poor,” the dimensions of the 
Hankel matrix should be “large” in order to get “good” error averaging. 

This brief exposure demonstrates that the Hankel approach imbeds a number of 
major drawbacks when addressing realization problems where the given set of Markov pa- 
rameters are contaminated by “significant” errors. 

On the same time, these drawbacks have stimulated additional research in an area 
indicated nowadays as “reduced order modelling” . The problem here is to approximate the 
given full-rank Hankel matrix by a low-rank matrix, which also has to be Hankel. 

For a number of applications [8], the latter stringent requirement might be relaxed to 
only finding a low-rank approximation, which not necessarily possesses the Hankel structure. 
For this approximate reduced-order modelling problem, a realization scheme also based on 
the QR factorization of the given Hankel matrix, is summarized in algorithm 4.1. Here the 
order n of the reduced-order system is assumed to be given. Furthermore, it is demonstrated 
that column pivoting might be incorporated in the QR factorization. 


Algorithm 4.1: 

STEP 1: Compute a QR factorization with column pivoting, indicated by the column per- 
mutation matrix n, of the given Hankel matrix Hi(M,N): 


H 1 (M,N)U^(Q i \Q 2 ) 


' Rn 

#12 \ 


0 

1*22 ) 



(27) 


From this decomposition, the Hankel matrix H i can be approximated by the matrix iff, 
given as: 

Hf = Q l \R u | Ru}n T = Q,[r cj | r c „ +) (28) 


with \\H i - Hf || 2 = HR 22 II 2 
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STEP 2: From the decomposition given in (28), the state space triplet {A^ , 6", (c“) T ) is 
then realized by solving the following least squares problems: 

min^«||Af [r c , ■ r CJV _ i ] - [r C2 • • • r e J|| 2 (29) 

min c » ||(c“) r [r Cl • • • r Cjv _ 1 ] - [fn •••h N _ 1 )|j 2 (30) 

and 6“ is given as: 

K = r ei (31) 

Equations (29-30) of algorithm 4.1 clearly demonstrate that the system triplet 
{A*,bc , c*} is realized via the solution of a set of least squares problems. This is com- 
pletely similar to the realization scheme based on the SVD [1]. On the same time, we 

observe from the same set of equations (29-30) that the condensed system structure of the 

realization is no longer preserved. Therefore, in comparison with algorithm 3.2, algorithm 
4.1 is no longer attractive , however it still remains elegant . This elegantness is due to the 
fact that algorithm 4.1 also only needs the upper triangular factor of the QR factorization. 

Furthermore, it is not necessary to incorporate column pivoting in algorithm 4.1. 
The necessity depends on the following trade-off. On the one hand, it makes the realization 
scheme computationally more complex , but on the other hand it allows the error in approx- 
imating the given Hankel matrix, as may be expressed by ||i? 22 ||, to be minimized. If, for 
example, the minimization of this error is crucial, it is recommended to use the column 
pivoting strategy described in [14] or [15]. 

In addition to the drawback of approximating the given Hankel matrix by a non- 
Hankel matrix (from which the system is realized), algorithm 4.1 also required the specifi- 
cation of the order of the reduced-order system. Since it was assumed that this information 
could not be retrieved from the second fundamental property of the Hankel approach, ad- 
ditional information is required. 

5. Experimental evaluation. In this section, algorithm 3.2 and 4.1 are 
evaluated experimentally. In this analysis, a third-order SISO system was taken as a test 
vehicle. The transfer function /?(^) of this plant is: 


13 



0.005496z 2 + 0.0202852 + 0.004672 
PK Z ) - z 3 _ 2.70066 z 2 + 2.424258Z - 0.72253 ^ ' 

and its impulse response is shown in Fig. 1. 

The experimental study begins with an evaluation of the size of the errors e t , for 
which algorithm 3.2 is still applicable. In judging the applicability we will evaluate the 
difference between the actual impulse response (calculated from (32)) and the one calculated 
from the realized state space triplet. 


5.1. Determining the applicability of algorithm 3.2. In this part of the 
evaluation, we furthermore compare algorithm 3.2 with the realization technique based on 
the SVD [1] and the one by Rissanen [3]. 


In a first experiment no additional errors on the Markov parameters, calculated 
from (32), are assumed. These Markov parameters were arranged in the Hankel matrix 
i?i(200,6). From this Hankel matrix a state space triplet {Aqr, bQR, Cq R } was realized 
via algorithm 3.2 and using the SVD as described in [l] the triplet { Asvd, bsvD, c ^vd) 
resulted. The impulse response calculated from both realizations coincides with the actual 
given one. This becomes clear from a comparison of the relative errors on the Markov 
parameters computed from each triplet. These relative errors are given as: 


INI» 


(33) 


and 


ll^tlb 


(34) 


for the triplet realized by the algorithm 3.2 and the one based on the SVD [l] respectively. 
Figure 2 graphically represents both these quantities and from this figure we clearly observe 
that the Markov parameters calculated from both realizations are of the same accuracy. 


In a second experiment, additional errors e t were considered. For the sequence {e*}, 
defined in (18), a zero-mean white-noise sequence with standard deviation a ei = 3 • 10 -6 is 
taken. This resulted in very small errors on the Markov parameters relative to their orig- 
inal magnitude. The Hankel matrix ffi(200,6) constructed from these perturbed Markov 
parameters is now used to realize two state space triplets based on the QR factorization 
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and the SVD. The relative errors (33) and (34) for this experimental condition are plotted 
in Fig. 3. From this figure, we observe that both relative errors are small and of the same 
order. However, we now observe from Fig. 3 that the results based on the SVD are more 
accurate. This discrepancy between the algorithm 3.2 and the method based on the SVD 
will increase when the size of the errors e x increases. Whether this discrepancy remains very 
small for the region where algorithm 3.2 remains applicable is investigated more closely in 
the next experiment. Let us first evaluate the performance of the technique developed by J. 
Rissanen [3] for the perturbed Markov parameters of the second experiment. The computed 
impulse response from this realization is plotted in combination with the one given in Fig. 4. 
This clearly reveals the severe limitations of this realization scheme. 

In the third experiment, we investigate the applicability of the order-detection mech- 
anism, summarized by equation (26). For this purpose, three different values of o ti were 
chosen, i.e. 10" 7 , 10" 6 and 10 -5 respectively. The information used in detecting whether 
the system is of order two is graphically represented in Fig. 5a. From this figure we clearly 
observe (the original system being third order) that a €t - 10" 5 determines the maximal size 
of the errors, for which (26) still detects the correct order. Of course the range of applica- 
bility of (26) could be somewhat extended for this application by decreasing the number of 
samples, however the size of the errors for which (26) detects the correct order, will remain 
relatively very small 

For a €i = 10“ 5 , the calculation of the impulse response based on the calculated 
triplet via algorithm 3.2 remains very small. The relative error (33) is again of the same 
order as those depicted in Fig. 3. This third experiment clarifies that the size of the errors 
e t for which algorithm 3.2 remains applicable have relatively to be very small. 

In Fig. 5b, the information necessary for determining whether the system is of order 
three is graphically represented. Here, obviously no restrictions imply. 

When the errors e, increase in magnitude, by making larer than 10" 5 , even the 
rank of the Hankel matrix becomes til-defined via the SVD. This is because no longer a 
“clear gap” is present in the sequence of the singular values of the “perturbed” Hankel 
matrix. Therefore, this confirms that when the size of the errors c x is outside the region 
of applicability of algorithm 3.2, also the scheme based on the SVD presented in [1] no 
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longer solves the partial realization problem straightforwardly. In this case, one focusses on 
the reduced- order modelling problem. The performance of algorithm 4.1 for this purpose is 
evaluated in the next section. 

5.2. Approximate reduced-order modelling. In this experiment we con- 
sidered errors e t outside the region of applicability of algorithm 3.2, and took <j C| — 10~ 3 . 
The detection of the rank based on finding a gap in the singular values of the constructed 
Hankel matrix tfi(100, 10) is not possible . Therefore, let us just postulate the order of 
the system. In this way, a system of order 3 and 9 were realized via algorithm 3.2. The 
computed impulse response from these realizations are respectively depicted in Fig. 6 and 
compared with the the given (noisy) impulse response. From Fig. 6 we clearly observe the 
bad performance of the third-order realized system and also the match of the ninth-order 
system is unsatisfactory. Another disadvantage of this last realization is although it is in 
its “condensed” controller Hessenberg form, its huge dimension (compared with the order 
of the actual system) will severely penalize the operation count of subsequent calculations 
based pn that realization. 

Using the same Hankel matrix Hi( 100, 10), we compared in, Fig. 7, the performance 
of algorithm 4.1 when no column pivoting was used with case where the column pivoting 
strategy described in [14] or [15] is used. In this same figure, we furthermore plotted the 
impulse response calculated from a realization obtained from a S VD of the given Hankel ma- 
trix. From this figure, we clearly observe the improvement and therefore also the necessity 
of using column pivoting in algorithm 4.1. 

6. The conclusions. In this paper, the use of the QR factorization is studied 
in solving the partial realization problem via the so-called Hankel approach [2]. 

Two different types of partial realization problems have been considered. 

In the first type, the Markov parameters are indicated to be “accurately” given, 
meaning that the errors on these parameters do not destroy the correspondence between 
the numerical rank of the Hankel matrix and the order of the system. 

For this type of realization problem, a new realization scheme has been proposed 
based on a straightforward QR factorization of the Hankel matrix. It has been clearly 
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demonstrated that this scheme outperforms the as “originally indicated computational at- 
trative schemes” developed by J. Rissanen [3], and L.S. De Jong [4]. This is because while 
the new scheme considers all the given Markov parameters and is due to the use of or- 
thogonal tranfomations numerically robust, it has also resulted in a realization which is 
in condensed system structure, i.e. the observer-Hessenberg or controller-Hessenberg form 
[9]. This allows the direct use of the vast amount of algorithms to solve in an efficient and 
reliable manner different system theoretical problems (see, for example, [9] for an overview). 

In an experimental comparison, it has been demonstrated that for the type of real- 
ization problem under consideration, the new scheme has a numerical reliability which is 
similar to the realization scheme based on the SVD [1]. However, the computational com- 
plexity of this realization scheme as well as the retrieval of a condensed system structure 
make this new scheme much more attractive. 

A second type of realization problems considered the errors on the Markov param- 
eters to violate the correspondence between the numerical rank of the Hankel matrix and 
the order of the system. For this type of realization problems, an algorithmic scheme was 
presented which approximated the given Hankel matrix by a low-rank non-Hankel matrix 
also based on the QR factorization. This type of realization problem was referred to in this 
paper as “approximate reduced-order modelling” in order to distinguish from the treated 
“reduced-order modelling problems” [5], [6], and [7], where the Hankel structure is imposed 
on the low-rank approximation. It was demonstrated here that column pivoting might 
be incorporated in this scheme. In comparison with the solution presented for the previ- 
ous type of realization problems, this second scheme (with column pivoting) possesses the 
same advantages over the SVD, except that now no condensed system structure is directly 
obtained. 

The main conclusion of this research is that the QR factorization with or without 
column pivoting might replace the SVD as computational tool in solving partial realization 
problems. The challenging question now is whether it might also help conceptually to better 
understand difficult areas, such as approximating the given Hankel matrix by a low-rank 
matrix, which also has the Hankel structure, in the partial realization problem. 
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Fig. 1: Impulse response of the third-order system with transfer function (32). 

Fig. 2: The relative error on the Markov parameters calculated from a state space triplet 
realized from a QR and SVD decomposition assuming noise-free Markov parameters given. 

Fig. 3: The relative error on the Markov parameters calculated from a state space triplet re- 
alized from a QR and SVD decomposition assuming perturbed Markov parameters 
(< 7 e , = 3-10 6 ). 

Fig. 4: Comparison of the impulse response realized via the technique developed by J. 
Rissanen with the given noisy impulse response (cr €| = 3 • 10 6 ). 

Fig. 5a: Evaluation of (26) in algorithm 3.2 in detecting whether the underlying system is 
of order two, for different error levels a €j . 

Fig. 5b: Evaluation of (26) in algorithm 3.2 in detecting whether the underlying system is 
of order three, for different error levels 

Fig. 6: The impulse response calculated from a state space model of orders 3 and 9 realized 
by algorithm 3.1. 

Fig. 7: The performance of algorithm 4.1 with and without column pivoting. 
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IMPULSE RESPONSE 


GIVEN IMPULSE RESPONSE 

USING ALGORITHM 3.2 (3 th ORDER SYSTEM) 

USING ALGORITHM 3.2 (9 th ORDER SYSTEM) 
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Fig. 6 


26 



IMPULSE RESPONSE 


GIVEN IMPULSE RESPONSE 

USING ALGORITHM 4.1 WITHOUT COLUMN PIVOTING 
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Fig. 7 
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